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RESUMEN 

Presentamos un modelo analitico de "capa delgada" de la interaccion de un 
chorro biconico y un medio ambiente en movimiento (alineado con la direction del 
chorro), asi como simulaciones numericas (axisimetricas) de dicha interaccion. Una 
situation similar, aunque en un escenario mas complejo, ocurre en la cabeza de la 
estructura cometaria de Mira. Por esta razon, en la mayoria de las simulaciones 
numericas exploramos parametros que son consistentes con las observaciones del 
flujo bipolar de Mira B. Para estos parametros la interaccion es no-radiativa, lo que 
resulta en una zona de interaccion jet /medio ambiente bast ante ancha. A pesar 
de esto, encontramos que el modelo analitico de capa delgada describe de manera 
satisfactoria la morfologia basica del flujo. 

ABSTRACT 

We present a "thin shell" model of the interaction of a biconical outflow and a 
streaming environment (aligned with the direction of the flow), as well as numerical 
(axisymmetric) simulations of such an interaction. A similar situation, although in 
a more complex setup, takes place at the head of the comet ary structure of Mira. 
Thus, for most of the numerical simulations we explore parameters consistent with 
the observed bipolar outflow from Mira B. For these parameters, the interaction is 
non-radiative, so that a rather broad jet/streaming environment interaction region 
is formed. In spite of this, a reasonable agreement between the thin-shell analytic 
model and the numerical simulations is obtained. 
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1. INTRODUCTION 

Prompted by the discovery of a striking ~ 2° 



- ISM: DYNAMICS — ISM: JETS 
STARS: AGB AND POST-AGB — 



Optical (Ha) observations by Meaburn et al. 



cometary tail attached to Mira by Mar tin et al. 
([20071), a series of theoretical and observational stud- 



ies of this object have been made. Theoretical works 



include 


Wareing et al. 


(2007); 


Raga et al. ( 


2008); 


Raga & Canto 


(2008); 


Esquivel et al. 


(2010), all of 



which are devoted to study the interaction between 
the post- AGB wind (from Mira A) and a streaming 
interstellar medium (ISM). This streaming medium 
corresponds to the motion of Mira through the galac- 
tic plane. Observations of the cometary tail have 
been done in several wavebands, including 21cm 
( |Matthews et al.||2008| ) and the IR ( |Ueta||2008 ). 
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(2009) revealed a fast bipolar outflow, which prob- 
ably arises from the less luminous component, Mira 
B. This outflow is embedded in the head of the 
cometary structure. It has radial velocities of 
±150 km s _1 (with respect to Mira), is relatively 
broad, and its projection onto the plane of the sky 
is approximately aligned with the direction of Mira's 
proper motion. 

The mass loss rate from Mira A has e stimated 
values from 10 -7 to 10 -5 M yr _1 (Knapp & Morris 
1985||Bowers fc Knapp|1988[ ). Part of this mass loss 
is accreted onto a disk around Mira B (and eventu- 
ally onto Mira B itself). Ireland et al. (2007) esti- 
mated a ~ 5 x 10 -9 M yr _1 accretion rate. Part of 
this accretion rate would then be redirected into the 
bipolar outflow from Mira B. 

In the Mira A/B system, we therefore appear to 
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have a bipolar outflow from Mira B which first in- 
teracts with the wind from Mira A and at larger 
distances emerges from the wind into the stream- 
ing ISM region. However, in the present paper we 
consider the more simple problem of an aligned (bi- 
conical) bipolar outflow/streaming environment in- 
teraction (i. e., neglecting the presence of the wind 
from Mira A). This is a first step to modeling the 
high velocity outflow observed by Meab um et al.| 
( [2009 ). 

We present an analytic, stationary, thin shell 
model of the interaction of a bi-conical outflow and 
an aligned, streaming environment in Section [2] 
Time-dependent, axisymmetric numerical simula- 
tions of this problem are described in Section [3J and 
their results presented in Section [4j With the sim- 
ulations we show the characteristics of the flow in 
greater detail, and evaluate the level of agreement 
with the analytic, thin shell solution. Finally, the 
implications of this model for the interpretation of 
Mira's outflow are discussed in Section |H 

2. THE ANALYTIC MODEL 

Let us consider a source that ejects a bipolar, con- 
ical outflow (with a half-opening angle 9j, velocity 
vj and mass loss rate Mj for each outflow lobe), im- 
mersed in a uniform, streaming environment of den- 
sity p a , which travels parallel to the outflow axis at a 
velocity v a . The situation is shown in the schematic 
diagram of Figure [TJ We assume that as the material 
from the outflow encounters the streaming environ- 
ment a thin shell of well mixed material is formed 
(the thick solid line shown in Fig[T]). This layer has 
a bow shock shape R(0), where R is the spherical 
radius and 6 the angle measured from the z-axis, see 
Figure 1). 
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Fig. 1. Schematic diagram showing the interaction be- 
tween a bi-conical outflow (of half-opening angle 0j , ve- 
locity Vj and total mass loss rate 2Mj) with a streaming 
environment of density p a and velocity v a (aligned with 
the outflow axis). 

At a position 6 along the bow shock, the material 
flowing along the thin shell has a mass rate: 

M{9) = / dMj + 7rr 2 p a v a , (1) 
Jo 



where the first term on the right represents the mass 
fed into the shell by the conical jet, and the second 
term is the mass from the streaming environment 
(with r = RsinQ being the cylindrical radius, see 
Figure [T]) . 

The radial momentum rate of the material in the 
shell is: 

n r (0) = M(O)v r (0) = [ vj sin dMj , (2) 
Jo 

with a contribution only from the jet. The axial 
momentum (along the z-direction) flowing along the 
shell is: 

r e 

U z (0) = M{0) v z (0) = / Vj cos dMj - irr 2 p a v\ . 
Jo 

(3) 

In equations ([2| and (J3|, v r and v z are the r- and 
z-components of the velocity of the (well mixed) ma- 
terial flowing along the thin shell. One can also write 
an expression for the angular momentum rate of the 
material flowing along the shell: 

J{0) = M(0)v e (0)R(0) = [ r'v a (2^p a v a )dr' 

Jo 

= ^7rp a v 2 a R 3 sin 3 0, (4) 

where vq is the velocity in the ^-direction (of the 
material flowing along the thin shell) , which is given 
by: 

vq = v r cos 6 — v z sin 6 . (5) 

Now, we consider a conical, bipolar outflow with 
a mass loss rate given by: 

f i M \ smOdO ; 6 < 6 j or > tt - j , 
dM j = <^ i -cos 04 - 3 ' 

[ ; j < < tt - j . 

(6) 

With this form for dMj, it is possible to carry out 
the integrals in equations ([l][3|, allowing to obtain 
the shape of the bow shock: 

R(6) = R cosec(9 x ^3 (1 - cotan 0) ; < 0j , 

(7) 

R(0) = R cosec 
x 3 [cos(# — 6j) sin 0j cosec — 0j cotan 6] ; 

0j < < tt - 0j , (8) 
R{6) = R cosec 

x ^3 [1 + (tt - - 20 j + sin 2!%) cotan 0) ; 

TT - 0j < < TT , (9) 



JET STREAMING ENVIRONMENT CAVITIES 



3 



where 



Ro = 



(10) 



27r(l - cos6j)p a vl ' 

is the on-axis standoff distance between the bow 
shock and the jet source. 




Fig. 2. Analytic solution for the shape of the bow shock 
produced by the interaction of a bi-conical outflow with 
an aligned, streaming environment. The outflow cones 



have a half- opening angle 



15°, and their outer 



boundaries are shown by the dashed lines. 

As an example, in Figure [2] we show the bow 
shock shape obtained from equations (Ffpl) for a bi- 



conical outflow with a half-opening angle 9j = 15°. 
The solution is quite fiat-topped in the < 0j "head" 
region. For > 0j the solution first rapidly curves, 
and then straightens out into a slowly expanding bow 
shock wing. 

3. THE NUMERICAL SIMULATIONS 

We have computed seven different axisymmet- 
ric gasdynamic simulations of aligned bi-conical 
jet /streaming environment interactions with the 
walicxe code, which is described in detail by |Es-| 
|quivel et al.| ( |2Q1Q ). The code integrates the gas- 
dynamic equations with an approximate Riemann 
solver (in this case we used a hybrid HLL-HLLC al- 
gorithm, see |Harten et al.|[l9g3l |Toro et al.| [19941 
Esquivel et al. 2010). A hydrogen ionization rate 



equation is solved along with the gasdynamic equa- 
tions in order to include the radiative losses through 
a parametrized cooling function, that depends on the 
density, temperature and hydrogen ionization frac- 
tion ( Raga fc Reipurth||2004 ). 

The walicxe code has a block based adaptive 
grid, with blocks of a constant number of cells (in 
our case, of 24 x 24 cells), which are refined by suc- 
cessive factors of 2. Most of our simulations were 
run at a medium resolution (labeled "mr" in Table 
1). These consist of two root blocks aligned in the 
axial direction, which are allowed to have 8 levels of 
refinement. The resolution at the finest level would 
correspond to 6144 x 3072 cells (axial and radial, re- 
spectively) in a uniform grid. In addition, we ran low 



and high resolution versions of model 1 (labeled "lr" 
and "hr" in Table 1). The low resolution run has 7 
levels of refinement (3072 x 1536 cells in a uniform 
grid), while the high resolution run has 9 levels of 
refinement (12288 x 6144 cells in a uniform grid). 

Models 1-6 have a computational box of the same 
size, [2 x 4]10 18 cm (radial and axial). Therefore 
the "mr" runs have a maximum resolution of 6.51 x 
10 14 cm, and the "lr" and "hr" runs have resolutions 
of 1.30 x 10 15 cm and 3.25 x 10 14 cm, respectively. 
In model 7 we were more interested in resolving the 
structure at the head of the shock than in the long 
cometary tail. To achieve a higher resolution (6.51 x 
10 13 cm) with the same number of computational 
cells, in these models we have reduced the size of the 
computational box by a factor of 10. 

The simulations were performed in a cylindrical 
computational grid with an axial extent [z min , z max ] 
(with the jet source placed at z = 0) and extend- 
ing radially from the symmetry axis (r = 0, where 
a reflecting boundary condition is applied) out to 
r max . An inflow condition is applied in the positive 
z grid boundary, where a streaming environment is 
injected with a constant velocity in the — z direc- 
tion. Outflow boundaries were applied at the neg- 
ative z grid boundary and at r max . In models 1-3 
Zmin = —2.8 x 10 18 cm and z max = 1.2 x 10 18 cm, 
in models 4-6 z m i n = —2.0 x 10 18 cm and z max = 
2.0 x 10 18 cm, in model 7 z min = —2.8 x 10 1T cm and 
z max = 1.2 x 10 17 cm. 

The parameters of the simulations are summa- 
rized in Table 1. The jets are imposed in two cones 
with a half-opening angle O and an initial outer ra- 
dius equal to Rj = 10 16 cm for models 1-6, and 
Rj = 10 15 cm for model 7. The outflows have a ve- 
locity Vj (see Table 1) and a temperature Tj = 100 K. 
The density inside the source follows a oc R~ 2 pro- 
file set by the mass loss rate Mj. The environment 
has an initial temperature T a — 10 4 K, a density 
n a and a velocity v a . The values for the parameters 
of models 1-6 are roughly consistent with the jet of 
Mira B while model 7 (with a denser streaming en- 
vironment) was chosen to produce a more radiative 
interaction region. Both the jet and the environ- 
ment are initially neutral, except for a seed electron 
density assumed to come from the presence of singly 
ionized carbon. 

In Table 1 we have also included the correspond- 
ing value of the bow-shock standoff distance Ro (see 
equation 10), and the cooling distances associated 
with the double shock structure at Ro. The cool- 
ing distance dj corresponds to the forward shock 
(with the streaming ambient medium), while d c r cor- 



4 



LOPEZ-CAMARA ET AL. 



responds to the cooling behind the reverse shock. In 
order to calculate these cooling distances, we have 
used the fit: 



100 cm -3 ' 



[8 x 10 



13 



cm I 



[3 x 10 n cm] 



100 km s" 



100 km s" 



-6.4 



(ii) 



to the cooling distances to 10 4 K given in the tabu- 
lation of self-consistent preionization, plane-parallel 



shock models of Hartigan et al. (1987). In equation 



(11), no is the pre-shock atom+ion number density, 



and uo is the shock velocity. 

4. RESULTS 

In Figures [3] -|5| we show the time-evolution (t = 
1 x 10 3 , 3xl0 3 , and 5xl0 3 yr, respectively) obtained 
from model 1 (see Table 1). For each time, we show 
the corresponding mesh configuration, the density 
nn (in cm -3 ) and the temperature T (in K) (top, 
middle, and bottom panels, respectively). Thanks 
to the adaptive grid of the code we used less than 
10% (for t x 10 3 yr) and 20% (for t > 10 3 yr) of the 
finest grid in the domain, thus reducing notably the 
required computer time. 

From the time-sequence shown in Figures [3] 
- [5} we see that the flow develops a cometary- 
shaped structure, due to the presence of the side- 
streaming environmental wind (which flows paral- 
lel to the z-axis, from right to left). The head of 
the cometary structure approaches the steady state 
star/bow shock stagnation distance Ro (see §2) at a 
time t « 10 3 yr. 

In Figure [6j we show single time frames (t = 5x 
10 3 yr) of the density and temperature stratifications 
(top and bottom parts, respectively, of each of the 
three panels) obtained from models 1 through 3 (see 
Table 1). In each plot we include the corresponding 
"thin shell" analytic solution (see Section [2]). 

Even though models 1 through 3 differ in Mj and 
6 (see Table 1), they result in basically the same 
flow morphology. Model 3 (which has the largest 
mass loss rate and opening angle) has the broad- 
est flow. Independently of the Mj and values, 
all three models coincide with the theoretical stag- 
nation point (Ro = 2.97 x 10 17 cm), and they also 
have tails with the same basic morphology. It is clear 
that the numerical simulations produce broad inter- 
action regions, which are a direct result of the fact 
that the post-shock cooling distances (d c r , and dj 
see Table 1) are not small compared to the width of 
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Fig. 3. Density (color scale given in cm -3 , middle panel) 
and temperature (color scale given in K, bottom panel) 
stratifications. The corresponding mesh configurations 
(on top of each density /temperature plot) are also shown, 
with each square representing a 24 x 24 grid point block. 
The results are for model 1 at t — 1 x 10 3 yr. The axes 
are labeled in units of 10 17 cm (the figure only shows the 
inner half of the radial extent, the rest is unperturbed.). 
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Fig. 4. Same as Figure [|] at t = 3 x 10 3 



yr. 
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Fig. 5. Same as Figure [|] at t = 5 x 10 3 yr. 
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TABLE 1 
MODEL CHARACTERISTICS. 



Model 


e 


V 3 


Mj 


V a 


n a 


Ro 


d c 7(i?otanB) 


d c f /(i? tanO) 


resolution 




[°] 


[km s- 1 ] 


[Mo s- 1 ] 


[km s- 1 ] 


[cm -3 ] 


[10 17 cm] 








1 


10 


200 


1.00 xlCT 10 


125 


0.044 


2.97 


402.4 


11.9 


lr,mr,hr 


2 


20 


200 


3.97xl0- 10 


125 


0.044 


2.97 


194.9 


5.7 


mr 


3 


30 


200 


8.82X10- 10 


125 


0.044 


2.97 


122.9 


3.6 


mr 


4 


10 


200 


1.00 xlCT 9 


125 


0.044 


9.39 


127.2 


3.7 


mr 


5 


20 


200 


3.97X10- 9 


125 


0.044 


9.39 


61.6 


1.8 


mr 


6 


30 


200 


8.82X10- 9 


125 


0.044 


9.39 


38.9 


1.1 


mr 


7 


10 


40 


l.OOxlO -10 


100 


1.0 


3.48 


0.28 


1.3 


mr 



the flows. We find that the thin layer analytic so- 
lution lies along a line dividing the computed flows 
into two regions: an internal low-density, hot region 
(corresponding to the jet cocoon) and a high-density 
external, cold region (corresponding to the shocked, 
streaming ambient medium). 

The fact that the cooling distances are large com- 
pared to the cross section of the jet indicates that the 
flows are non-radiative. In this case it is remarkable 
that the analytical solution does successfully divide 
the hot /low-density from the cold/high-density re- 
gions. This is not necessarily expected, as the ana- 
lytical model assumes that the material is well mixed 
in a thin layer, which would arise more naturally in 
radiative flows (where cooling occurs in a small re- 
gion). 

The density and temperature stratifications (at 
t = 5 x 10 3 yr) for models 4-6 are shown in Fig- 
ure [7] These models, with mass loss rates one order 
of magnitude larger than those from the first three 
models (~ 10 -9 M s _1 , for details see Table 1) have 
smaller cooling distances than models 1-3, however 
they are still in the non-radiative regime. The three 
models roughly match the analytically derived stag- 
nation point (Ro = 9.39 x 10 17 cm), and produce 
broader flows for larger values of the opening angle 
6. 

In order to check how well the numerical simu- 
lations reproduced the theoretical stagnation point 
in the thin shell analytical solution from Section [2j 
we show a zoom of the density structures of the jet 
heads for all models in Figures [8][9j We find that 
a steady state is reached (with small discrepancies, 
more evident in models 4-7, Fig. |9|, with a shock 
structure approximately centered on the stagnation 
radius Ro (see Table 1). One can also see that the 
place at which the solution curves sharply coincides 
nicely with the end of the dense wall of the jet. 

In order to illustrate the effect of the numerical 
resolution, we computed two extra simulations with 
the setup of model 1. The resulting density strati- 
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Fig. 6. Density stratifications (color scale given in cm -3 , 
blue panels); and temperature (color scale given in K, or- 
ange panels). Top, middle and bottom figures correspond 
to models 1, 2 and 3 (respectively). The integration time 
is t=5xl0 3 yrs (top). The axes are labeled in units of 



fications for an integration time t = 5 x 10 3 yr, are 



shown in Figure 10 



Even though more complex structures are ob- 
tained for increasing resolutions, the basic morphol- 
ogy of the flow (stagnation point, cocoon size, den- 
sity profile, etc.) is the same in all cases. Taking the 
number of grid points across the jet as an estimate 
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Fig. 7. Same as Figure [6] but for models 4, 5, and 6 (top, 
middle, and bottom respectively). 

of the Reynolds number of the simulation (Re), in 
these models we only reach Re ~ Rj/Ax « 10 2 (in 
our hr run). To avoid Reynolds number dependency, 
the simulations should have resolutions larger by a 
factor of at least two orders of magnitude greater, 
which is presently unattainable. However, from the 
results shown in Figure 8, it is clear that the mr 
resolution appears to produce the correct basic mor- 
phology for the flow. We must note that this re- 
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suit is similar to the one found by Wareing et al 



(2007). In this, the authors simulated the interaction 



Fig. 8. Density stratifications (color scale given in cm -3 ). 
Top, middle and bottom figures correspond to model 1 
(and 2); 3; and 4 (respectively). The integration time 
is t=5xl0 3 yrs (top). The axes are labeled in units of 
10 17 cm. The broken line delimits the analytic solution 
proposed in Section [2] 

of an AGB star moving through the ISM (with also 
Re - 10 2 ). The key result was that independently of 
the Reynolds number, the ISM ram pressure stripped 
material from the head of the AGBOs bow shock. 
Such stripping as well as the presence of Rayleigh- 
Taylor instabilities produced turbulent clump struc- 
tures that eventually moved down stream to the tail 
of the bow shock. This mechanism even though in 
smaller scale, is also present in our simulations (see 
Figures |3]-[5|. 

As we have mentioned earlier, the analytic model 
assumes that the material of the environment and 
the jet mix well in a thin shell. Such a situation 
is more easily met in radiative flows, characterized 
by a thin cooling layer. However, our simulations 
with parameters appropriate for the jet of Mir a B 
are basically non-radiative. In order to check what 
happens with a more radiative flow we ran model 
7. In this model, the reverse shock is indeed radia- 
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and therefore a 10 times larger resolution than the 
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Fig. 9. Same as Figure [8] but for models 5; 6; and 7 
middle and bottom respectively). 
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Fig. 10. Density stratifications (color scale given in 
cm -3 ) for model 1 at t=3xl0 3 yrs, for an "hr" reso- 
lution (top panel); "mr" resolution (middle panel); and 
"lr" resolution (bottom panel). The axes are labeled in 
units of 10 17 cm. 

tive (d c r < Rq tan 0) while the cooling distance of 
the forward shock is of the order of the cross sec- 
tion of the jet (dj ~ i?otanO) (see Table 1). This 
model has a 10 times smaller computational domain, 



mr models (see Section 3). In Figure 11, we overlap 
the thin layer analytic solution on top of the den- 
sity and temperature stratifications obtained from 
model 7 at various times (v a = 100 km s _1 , and 
Vj = 40 km s _1 ). For this model the stagnation point 
is Ro = 3.48 x 10 16 cm, and the cooling distances are 
d c r = 0.28 Ro tan 9, and dj = 1.31i? tan9. As 
can be seen from the time-sequence, a steady state 
is not reached in this model. The model develops a 
turbulent structure which is not stable, and the dis- 
tance from the jet source to the bow-shock oscillates 
around Rq. 
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Fig. 11. Density stratifications (color scale given in 
cm -3 , blue panels); and temperature (color scale given 
in K, orange panels) of model 7 at t=10xl0 3 yrs (top); 
t=18xl0 3 yrs (middle); and t=24xl0 3 yrs (bottom). 
The axes are labeled in units of 10 16 cm. 

5. DISCUSION 

Motivated by the recently observed bipolar jet 
from Mira ( Meaburn et al.|2009 ), we study the inter- 
action of a steady bi-conical outflow with a stream- 
ing environment. We use this scenario to perform 2D 
axisymmetrical simulations and an analytical model 
of such an interaction. 
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We have studied the particular case in which the 
bi-conical jets are aligned with the streaming envi- 
ronment (as shown in Figure [TJ. This of course is 
an idealization of the situation which could possibly 
be found in Mira, which is unlikely to have a perfect 
alignment between the outflow axis and the direc- 
tion of Mira's motion. It allows, however, simple 
axisymmetrical analytical and numerical solutions 
which are useful as a first exploration of the more 
complex, less symmetric problem of the jet in Mira 
B. 

The proposed "thin shell" analytical solution, 
based on the momentum equilibrium between the 
outflow material from the conical jet and the am- 
bient medium, traces the shape of the shock (see 
Equations 7-9) . The model assumes that the interac- 
tion occurs in a very thin layer in which the material 
from the jet and the ISM mixes well. This situation 
is likely to occur in radiative flows, there the size of 
the cooling distances are small compared with the 
flow dimensions. This, however is not the case of the 
resulting interaction of the outflow launched from 
Mira B and its surrounding environment. 

The models reach a stationary state with a fixed 
stand-off distance between the jet head and the out- 
flow source. This configuration might be relevant for 
the collimated outflow from Mira, in which the for- 
ward directed jet head has a very low proper motion 



the binary Mira AB system, the mass loss rate from 
Mira A's uncollimat ed wind (with M ~10~ 6 M©s _1 , 
Ireland et al. (2007)) clearly dominates the interac- 



with respect to Mira (see |Meaburn et aL (2009)). 



We find that the non-radiative flow configura- 
tions (that result from the parameters deduced from 
the bipolar outflow from Mira) have broad interac- 
tion regions which approximately coincide with the 
predictions from the analytic model (in which a thin 
shell interaction is assumed). Also, the analytical 
solutions lie close to the contact discontinuity divid- 
ing the material in the jet cocoon from the shocked 
environment region. 

We have tested the convergence of the numeri- 
cal simulations for increasing resolutions. We find 
that at different resolutions we obtain similar large- 
scale flow structures, but with different structures 
at smaller scales (in which the Kelvin-Helmholtz in- 
stabilities associated with shear layers are active). 
Better convergence is expected only for resolutions 
higher by at least two orders of magnitude, in which 
the simulations would start to approach the "high 
Reynolds number regime" appropriate for the real, 
astrophysical flow. 

We end by noting again that the present model of 
a bipolar jet/streaming ISM interaction is not com- 
pletely consistent with the outflow from the Mira 
system ( |Meaburn et al.| 2009). This is because in 



tion with the ambient medium. The mass loss from 
the bi-conical jet (Mj ~lO _9 M s _1 , inferred from 



Martin et al. ( 2007[ )), will simply modify the mor- 
phology created by the interaction from the isotropic 
wind and the ambient medium, except in the region 
in which the leading jet from Mira B emerges from 
the Mira A wind bubble. This is the region to which 
our model would in principle apply. 

Future models of the outflow from Mira should 
include both the wind from the primary star and the 
collimated outflow from the secondary. This more 
complex flow will clearly be the topic of future papers 
on Mira's outflow/ISM interaction. 
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